{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Source of the materials**: Biopython cookbook (adapted)\n",
    "<font color='red'>Status: Draft</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BLAST\n",
    "===\n",
    "[Running BLAST over the Internet](#Running-BLAST-over-the-Internet)\n",
    "\n",
    "[Saving blast output](#Saving-blast-output)\n",
    "\n",
    "[Running BLAST locally](#Running-BLAST-locally)\n",
    "\n",
    "[Parsing BLAST output](#Parsing-BLAST-output)\n",
    "\n",
    "[The BLAST record class](#The-BLAST-record-class)\n",
    "\n",
    "[Parsing plain-text BLAST output](#Deprecated-BLAST-parsers)\n",
    "\n",
    "Hey, everybody loves BLAST right? I mean, geez, how can it get any\n",
    "easier to do comparisons between one of your sequences and every other\n",
    "sequence in the known world? But, of course, this section isn’t about\n",
    "how cool BLAST is, since we already know that. It is about the problem\n",
    "with BLAST – it can be really difficult to deal with the volume of data\n",
    "generated by large runs, and to automate BLAST runs in general.\n",
    "\n",
    "Fortunately, the Biopython folks know this only too well, so they’ve\n",
    "developed lots of tools for dealing with BLAST and making things much\n",
    "easier. This section details how to use these tools and do useful things\n",
    "with them.\n",
    "\n",
    "Dealing with BLAST can be split up into two steps, both of which can be\n",
    "done from within Biopython. Firstly, running BLAST for your query\n",
    "sequence(s), and getting some output. Secondly, parsing the BLAST output\n",
    "in Python for further analysis.\n",
    "\n",
    "Your first introduction to running BLAST was probably via the NCBI\n",
    "web-service. In fact, there are lots of ways you can run BLAST, which\n",
    "can be categorised in several ways. The most important distinction is\n",
    "running BLAST locally (on your own machine), and running BLAST remotely\n",
    "(on another machine, typically the NCBI servers). We’re going to start\n",
    "this chapter by invoking the NCBI online BLAST service from within a\n",
    "Python script.\n",
    "\n",
    "*NOTE*: The following Chapter \\[chapter:searchio\\] describes\n",
    "`Bio.SearchIO`, an *experimental* module in Biopython. We intend this to\n",
    "ultimately replace the older `Bio.Blast` module, as it provides a more\n",
    "general framework handling other related sequence searching tools as\n",
    "well. However, until that is declared stable, for production code please\n",
    "continue to use the `Bio.Blast` module for dealing with NCBI BLAST."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running BLAST over the Internet\n",
    "-----------------------------\n",
    "\n",
    "We use the function `qblast()` in the `Bio.Blast.NCBIWWW` module to call\n",
    "the online version of BLAST. This has three non-optional arguments:\n",
    "\n",
    "-   The first argument is the blast program to use for the search, as a\n",
    "    lower case string. The options and descriptions of the programs are\n",
    "    available at\n",
    "    <http://www.ncbi.nlm.nih.gov/BLAST/blast_program.shtml>. Currently\n",
    "    `qblast` only works with blastn, blastp, blastx, tblast and tblastx.\n",
    "\n",
    "-   The second argument specifies the databases to search against.\n",
    "    Again, the options for this are available on the NCBI web pages at\n",
    "    <http://www.ncbi.nlm.nih.gov/BLAST/blast_databases.shtml>.\n",
    "\n",
    "-   The third argument is a string containing your query sequence. This\n",
    "    can either be the sequence itself, the sequence in fasta format, or\n",
    "    an identifier like a GI number.\n",
    "\n",
    "The `qblast` function also take a number of other option arguments which\n",
    "are basically analogous to the different parameters you can set on the\n",
    "BLAST web page. We’ll just highlight a few of them here:\n",
    "\n",
    "-   The `qblast` function can return the BLAST results in various\n",
    "    formats, which you can choose with the optional `format_type`\n",
    "    keyword: `\"HTML\"`, `\"Text\"`, `\"ASN.1\"`, or `\"XML\"`. The default is\n",
    "    `\"XML\"`, as that is the format expected by the parser, described in\n",
    "    section \\[sec:parsing-blast\\] below.\n",
    "\n",
    "-   The argument `expect` sets the expectation or e-value threshold.\n",
    "\n",
    "For more about the optional BLAST arguments, we refer you to the NCBI’s\n",
    "own documentation, or that built into Biopython:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function qblast in module Bio.Blast.NCBIWWW:\n",
      "\n",
      "qblast(program, database, sequence, auto_format=None, composition_based_statistics=None, db_genetic_code=None, endpoints=None, entrez_query='(none)', expect=10.0, filter=None, gapcosts=None, genetic_code=None, hitlist_size=50, i_thresh=None, layout=None, lcase_mask=None, matrix_name=None, nucl_penalty=None, nucl_reward=None, other_advanced=None, perc_ident=None, phi_pattern=None, query_file=None, query_believe_defline=None, query_from=None, query_to=None, searchsp_eff=None, service=None, threshold=None, ungapped_alignment=None, word_size=None, alignments=500, alignment_view=None, descriptions=500, entrez_links_new_window=None, expect_low=None, expect_high=None, format_entrez_query=None, format_object=None, format_type='XML', ncbi_gi=None, results_file=None, show_overview=None, megablast=None)\n",
      "    Do a BLAST search using the QBLAST server at NCBI.\n",
      "    \n",
      "    Supports all parameters of the qblast API for Put and Get.\n",
      "    Some useful parameters:\n",
      "    \n",
      "     - program        blastn, blastp, blastx, tblastn, or tblastx (lower case)\n",
      "     - database       Which database to search against (e.g. \"nr\").\n",
      "     - sequence       The sequence to search.\n",
      "     - ncbi_gi        TRUE/FALSE whether to give 'gi' identifier.\n",
      "     - descriptions   Number of descriptions to show.  Def 500.\n",
      "     - alignments     Number of alignments to show.  Def 500.\n",
      "     - expect         An expect value cutoff.  Def 10.0.\n",
      "     - matrix_name    Specify an alt. matrix (PAM30, PAM70, BLOSUM80, BLOSUM45).\n",
      "     - filter         \"none\" turns off filtering.  Default no filtering\n",
      "     - format_type    \"HTML\", \"Text\", \"ASN.1\", or \"XML\".  Def. \"XML\".\n",
      "     - entrez_query   Entrez query to limit Blast search\n",
      "     - hitlist_size   Number of hits to return. Default 50\n",
      "     - megablast      TRUE/FALSE whether to use MEga BLAST algorithm (blastn only)\n",
      "     - service        plain, psi, phi, rpsblast, megablast (lower case)\n",
      "    \n",
      "    This function does no checking of the validity of the parameters\n",
      "    and passes the values to the server as is.  More help is available at:\n",
      "    http://www.ncbi.nlm.nih.gov/BLAST/Doc/urlapi.html\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from Bio.Blast import NCBIWWW\n",
    "help(NCBIWWW.qblast)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Note that the default settings on the NCBI BLAST website are not quite\n",
    "the same as the defaults on QBLAST. If you get different results, you’ll\n",
    "need to check the parameters (e.g., the expectation value threshold and\n",
    "the gap values).\n",
    "\n",
    "For example, if you have a nucleotide sequence you want to search\n",
    "against the nucleotide database (nt) using BLASTN, and you know the GI\n",
    "number of your query sequence, you can use:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "URLError",
     "evalue": "<urlopen error [Errno -3] Temporary failure in name resolution>",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mgaierror\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py\u001b[0m in \u001b[0;36mdo_open\u001b[1;34m(self, http_class, req, **http_conn_args)\u001b[0m\n\u001b[0;32m   1239\u001b[0m             \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1240\u001b[1;33m                 \u001b[0mh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrequest\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mreq\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_method\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mreq\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mselector\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mreq\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mheaders\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1241\u001b[0m             \u001b[1;32mexcept\u001b[0m \u001b[0mOSError\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0merr\u001b[0m\u001b[1;33m:\u001b[0m \u001b[1;31m# timeout error\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/http/client.py\u001b[0m in \u001b[0;36mrequest\u001b[1;34m(self, method, url, body, headers)\u001b[0m\n\u001b[0;32m   1082\u001b[0m         \u001b[1;34m\"\"\"Send a complete request to the server.\"\"\"\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1083\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_send_request\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmethod\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0murl\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbody\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mheaders\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1084\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/http/client.py\u001b[0m in \u001b[0;36m_send_request\u001b[1;34m(self, method, url, body, headers)\u001b[0m\n\u001b[0;32m   1127\u001b[0m             \u001b[0mbody\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mbody\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mencode\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'iso-8859-1'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1128\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mendheaders\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbody\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1129\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/http/client.py\u001b[0m in \u001b[0;36mendheaders\u001b[1;34m(self, message_body)\u001b[0m\n\u001b[0;32m   1078\u001b[0m             \u001b[1;32mraise\u001b[0m \u001b[0mCannotSendHeader\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1079\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_send_output\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmessage_body\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1080\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/http/client.py\u001b[0m in \u001b[0;36m_send_output\u001b[1;34m(self, message_body)\u001b[0m\n\u001b[0;32m    910\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 911\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    912\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mmessage_body\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/http/client.py\u001b[0m in \u001b[0;36msend\u001b[1;34m(self, data)\u001b[0m\n\u001b[0;32m    853\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mauto_open\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 854\u001b[1;33m                 \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mconnect\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    855\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/http/client.py\u001b[0m in \u001b[0;36mconnect\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    825\u001b[0m         self.sock = self._create_connection(\n\u001b[1;32m--> 826\u001b[1;33m             (self.host,self.port), self.timeout, self.source_address)\n\u001b[0m\u001b[0;32m    827\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msock\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msetsockopt\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msocket\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mIPPROTO_TCP\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msocket\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mTCP_NODELAY\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/socket.py\u001b[0m in \u001b[0;36mcreate_connection\u001b[1;34m(address, timeout, source_address)\u001b[0m\n\u001b[0;32m    692\u001b[0m     \u001b[0merr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 693\u001b[1;33m     \u001b[1;32mfor\u001b[0m \u001b[0mres\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mgetaddrinfo\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mhost\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mport\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mSOCK_STREAM\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    694\u001b[0m         \u001b[0maf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msocktype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mproto\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcanonname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msa\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/socket.py\u001b[0m in \u001b[0;36mgetaddrinfo\u001b[1;34m(host, port, family, type, proto, flags)\u001b[0m\n\u001b[0;32m    731\u001b[0m     \u001b[0maddrlist\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 732\u001b[1;33m     \u001b[1;32mfor\u001b[0m \u001b[0mres\u001b[0m \u001b[1;32min\u001b[0m \u001b[0m_socket\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetaddrinfo\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mhost\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mport\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfamily\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mproto\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    733\u001b[0m         \u001b[0maf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msocktype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mproto\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcanonname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msa\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mgaierror\u001b[0m: [Errno -3] Temporary failure in name resolution",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[1;31mURLError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-2-f79554f7040a>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mBio\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mBlast\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mNCBIWWW\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mresult_handle\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mNCBIWWW\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mqblast\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"blastn\"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"nt\"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"8332116\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/site-packages/Bio/Blast/NCBIWWW.py\u001b[0m in \u001b[0;36mqblast\u001b[1;34m(program, database, sequence, auto_format, composition_based_statistics, db_genetic_code, endpoints, entrez_query, expect, filter, gapcosts, genetic_code, hitlist_size, i_thresh, layout, lcase_mask, matrix_name, nucl_penalty, nucl_reward, other_advanced, perc_ident, phi_pattern, query_file, query_believe_defline, query_from, query_to, searchsp_eff, service, threshold, ungapped_alignment, word_size, alignments, alignment_view, descriptions, entrez_links_new_window, expect_low, expect_high, format_entrez_query, format_object, format_type, ncbi_gi, results_file, show_overview, megablast)\u001b[0m\n\u001b[0;32m    164\u001b[0m                            \u001b[0mmessage\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    165\u001b[0m                            {\"User-Agent\": \"BiopythonClient\"})\n\u001b[1;32m--> 166\u001b[1;33m         \u001b[0mhandle\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_urlopen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrequest\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    167\u001b[0m         \u001b[0mresults\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_as_string\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mhandle\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mread\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    168\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py\u001b[0m in \u001b[0;36murlopen\u001b[1;34m(url, data, timeout, cafile, capath, cadefault, context)\u001b[0m\n\u001b[0;32m    160\u001b[0m     \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    161\u001b[0m         \u001b[0mopener\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_opener\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 162\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mopener\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0murl\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    163\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    164\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0minstall_opener\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mopener\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py\u001b[0m in \u001b[0;36mopen\u001b[1;34m(self, fullurl, data, timeout)\u001b[0m\n\u001b[0;32m    463\u001b[0m             \u001b[0mreq\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmeth\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mreq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    464\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 465\u001b[1;33m         \u001b[0mresponse\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_open\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mreq\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    466\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    467\u001b[0m         \u001b[1;31m# post-process response\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py\u001b[0m in \u001b[0;36m_open\u001b[1;34m(self, req, data)\u001b[0m\n\u001b[0;32m    481\u001b[0m         \u001b[0mprotocol\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mreq\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtype\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    482\u001b[0m         result = self._call_chain(self.handle_open, protocol, protocol +\n\u001b[1;32m--> 483\u001b[1;33m                                   '_open', req)\n\u001b[0m\u001b[0;32m    484\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    485\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py\u001b[0m in \u001b[0;36m_call_chain\u001b[1;34m(self, chain, kind, meth_name, *args)\u001b[0m\n\u001b[0;32m    441\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mhandler\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mhandlers\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    442\u001b[0m             \u001b[0mfunc\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mhandler\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmeth_name\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 443\u001b[1;33m             \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    444\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mresult\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    445\u001b[0m                 \u001b[1;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py\u001b[0m in \u001b[0;36mhttp_open\u001b[1;34m(self, req)\u001b[0m\n\u001b[0;32m   1266\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1267\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mhttp_open\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mreq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1268\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdo_open\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mhttp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mclient\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mHTTPConnection\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mreq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1269\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1270\u001b[0m     \u001b[0mhttp_request\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mAbstractHTTPHandler\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdo_request_\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py\u001b[0m in \u001b[0;36mdo_open\u001b[1;34m(self, http_class, req, **http_conn_args)\u001b[0m\n\u001b[0;32m   1240\u001b[0m                 \u001b[0mh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrequest\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mreq\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_method\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mreq\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mselector\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mreq\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mheaders\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1241\u001b[0m             \u001b[1;32mexcept\u001b[0m \u001b[0mOSError\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0merr\u001b[0m\u001b[1;33m:\u001b[0m \u001b[1;31m# timeout error\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1242\u001b[1;33m                 \u001b[1;32mraise\u001b[0m \u001b[0mURLError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1243\u001b[0m             \u001b[0mr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetresponse\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1244\u001b[0m         \u001b[1;32mexcept\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mURLError\u001b[0m: <urlopen error [Errno -3] Temporary failure in name resolution>"
     ]
    }
   ],
   "source": [
    "from Bio.Blast import NCBIWWW\n",
    "result_handle = NCBIWWW.qblast(\"blastn\", \"nt\", \"8332116\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Alternatively, if we have our query sequence already in a FASTA\n",
    "formatted file, we just need to open the file and read in this record as\n",
    "a string, and use that as the query argument:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIWWW\n",
    "fasta_string = open(\"data/m_cold.fasta\").read()\n",
    "result_handle = NCBIWWW.qblast(\"blastn\", \"nt\", fasta_string)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We could also have read in the FASTA file as a `SeqRecord` and then\n",
    "supplied just the sequence itself:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIWWW\n",
    "from Bio import SeqIO\n",
    "record = SeqIO.read(\"data/m_cold.fasta\", format=\"fasta\")\n",
    "result_handle = NCBIWWW.qblast(\"blastn\", \"nt\", record.seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Supplying just the sequence means that BLAST will assign an identifier\n",
    "for your sequence automatically. You might prefer to use the `SeqRecord`\n",
    "object’s format method to make a FASTA string (which will include the\n",
    "existing identifier):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIWWW\n",
    "from Bio import SeqIO\n",
    "record = SeqIO.read(\"data/m_cold.fasta\", format=\"fasta\")\n",
    "result_handle = NCBIWWW.qblast(\"blastn\", \"nt\", record.format(\"fasta\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This approach makes more sense if you have your sequence(s) in a\n",
    "non-FASTA file format which you can extract using `Bio.SeqIO` (see\n",
    "Chapter [5 - Sequence Input and Output](05 - Sequence Input and Output.ipynb).)\n",
    "\n",
    "Whatever arguments you give the `qblast()` function, you should get back\n",
    "your results in a handle object (by default in XML format). The next\n",
    "step would be to parse the XML output into Python objects representing\n",
    "the search results (Section \\[sec:parsing-blast\\]), but you might want\n",
    "to save a local copy of the output file first. I find this especially\n",
    "useful when debugging my code that extracts info from the BLAST results\n",
    "(because re-running the online search is slow and wastes the NCBI\n",
    "computer time)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Saving blast output\n",
    "------------------\n",
    "\n",
    "We need to be a bit careful since we can use `result_handle.read()` to\n",
    "read the BLAST output only once – calling `result_handle.read()` again\n",
    "returns an empty string."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with open(\"my_blast.xml\", \"w\") as save_to:\n",
    "    save_to.write(result_handle.read())\n",
    "    result_handle.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "After doing this, the results are in the file `my_blast.xml` and the\n",
    "original handle has had all its data extracted (so we closed it).\n",
    "However, the `parse` function of the BLAST parser (described\n",
    "in \\[sec:parsing-blast\\]) takes a file-handle-like object, so we can\n",
    "just open the saved file for input:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with open(\"my_blast.xml\") as result_handle:\n",
    "    print(result_handle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we’ve got the BLAST results back into a handle again, we are\n",
    "ready to do something with them, so this leads us right into the parsing\n",
    "section (see Section \\[sec:parsing-blast\\] below). You may want to jump\n",
    "ahead to that now …."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running BLAST locally\n",
    "--------------------\n",
    "\n",
    "### Introduction\n",
    "\n",
    "Running BLAST locally (as opposed to over the internet, see\n",
    "Section \\[sec:running-www-blast\\]) has at least major two advantages:\n",
    "\n",
    "-   Local BLAST may be faster than BLAST over the internet;\n",
    "\n",
    "-   Local BLAST allows you to make your own database to search for\n",
    "    sequences against.\n",
    "\n",
    "Dealing with proprietary or unpublished sequence data can be another\n",
    "reason to run BLAST locally. You may not be allowed to redistribute the\n",
    "sequences, so submitting them to the NCBI as a BLAST query would not be\n",
    "an option.\n",
    "\n",
    "Unfortunately, there are some major drawbacks too – installing all the\n",
    "bits and getting it setup right takes some effort:\n",
    "\n",
    "-   Local BLAST requires command line tools to be installed.\n",
    "\n",
    "-   Local BLAST requires (large) BLAST databases to be setup (and\n",
    "    potentially kept up to date).\n",
    "\n",
    "To further confuse matters there are several different BLAST packages\n",
    "available, and there are also other tools which can produce imitation\n",
    "BLAST output files, such as BLAT.\n",
    "\n",
    "### Standalone NCBI BLAST+\n",
    "\n",
    "The “new” [NCBI\n",
    "BLAST+](http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=Download)\n",
    "suite was released in 2009. This replaces the old NCBI “legacy” BLAST\n",
    "package (see below).\n",
    "\n",
    "This section will show briefly how to use these tools from within\n",
    "Python. If you have already read or tried the alignment tool examples in\n",
    "Section \\[sec:alignment-tools\\] this should all seem quite\n",
    "straightforward. First, we construct a command line string (as you would\n",
    "type in at the command line prompt if running standalone BLAST by hand).\n",
    "Then we can execute this command from within Python.\n",
    "\n",
    "For example, taking a FASTA file of gene nucleotide sequences, you might\n",
    "want to run a BLASTX (translation) search against the non-redundant (NR)\n",
    "protein database. Assuming you (or your systems administrator) has\n",
    "downloaded and installed the NR database, you might run:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "blastx -query opuntia.fasta -db nr -out opuntia.xml -evalue 0.001 -outfmt 5\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This should run BLASTX against the NR database, using an expectation\n",
    "cut-off value of $0.001$ and produce XML output to the specified file\n",
    "(which we can then parse). On my computer this takes about six minutes -\n",
    "a good reason to save the output to a file so you can repeat any\n",
    "analysis as needed.\n",
    "\n",
    "From within Biopython we can use the NCBI BLASTX wrapper from the\n",
    "`Bio.Blast.Applications` module to build the command line string, and\n",
    "run it:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from Bio.Blast.Applications import NcbiblastxCommandline\n",
    "help(NcbiblastxCommandline)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "blastx_cline = NcbiblastxCommandline(query=\"data/opuntia.fasta\", db=\"nr\", evalue=0.001,\n",
    "outfmt=5, out=\"opuntia.xml\")\n",
    "blastx_cline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(blastx_cline)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>TODO: Need to add protein database [nr] or change example</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# stdout, stderr = blastx_cline()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example there shouldn’t be any output from BLASTX to the\n",
    "terminal, so stdout and stderr should be empty. You may want to check\n",
    "the output file `opuntia.xml` has been created.\n",
    "\n",
    "As you may recall from earlier examples in the tutorial, the\n",
    "`opuntia.fasta` contains seven sequences, so the BLAST XML output should\n",
    "contain multiple results. Therefore use `Bio.Blast.NCBIXML.parse()` to\n",
    "parse it as described below in Section \\[sec:parsing-blast\\].\n",
    "\n",
    "### Other versions of BLAST\n",
    "\n",
    "NCBI BLAST+ (written in C++) was first released in 2009 as a replacement\n",
    "for the original NCBI “legacy” BLAST (written in C) which is no longer\n",
    "being updated. There were a lot of changes – the old version had a\n",
    "single core command line tool `blastall` which covered multiple\n",
    "different BLAST search types (which are now separate commands in\n",
    "BLAST+), and all the command line options were renamed. Biopython’s\n",
    "wrappers for the NCBI “legacy” BLAST tools have been deprecated and will\n",
    "be removed in a future release. To try to avoid confusion, we do not\n",
    "cover calling these old tools from Biopython in this tutorial.\n",
    "\n",
    "You may also come across [Washington University\n",
    "BLAST](http://blast.wustl.edu/) (WU-BLAST), and its successor, [Advanced\n",
    "Biocomputing BLAST](http://blast.advbiocomp.com) (AB-BLAST, released in\n",
    "2009, not free/open source). These packages include the command line\n",
    "tools `wu-blastall` and `ab-blastall`, which mimicked `blastall` from\n",
    "the NCBI “legacy” BLAST suite. Biopython does not currently provide\n",
    "wrappers for calling these tools, but should be able to parse any NCBI\n",
    "compatible output from them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parsing BLAST output\n",
    "-------------------\n",
    "\n",
    "As mentioned above, BLAST can generate output in various formats, such\n",
    "as XML, HTML, and plain text. Originally, Biopython had parsers for\n",
    "BLAST plain text and HTML output, as these were the only output formats\n",
    "offered at the time. Unfortunately, the BLAST output in these formats\n",
    "kept changing, each time breaking the Biopython parsers. Our HTML BLAST\n",
    "parser has been removed, but the plain text BLAST parser is still\n",
    "available (see Section \\[sec:parsing-blast-deprecated\\]). Use it at your\n",
    "own risk, it may or may not work, depending on which BLAST version\n",
    "you’re using.\n",
    "\n",
    "As keeping up with changes in BLAST became a hopeless endeavor,\n",
    "especially with users running different BLAST versions, we now recommend\n",
    "to parse the output in XML format, which can be generated by recent\n",
    "versions of BLAST. Not only is the XML output more stable than the plain\n",
    "text and HTML output, it is also much easier to parse automatically,\n",
    "making Biopython a whole lot more stable.\n",
    "\n",
    "You can get BLAST output in XML format in various ways. For the parser,\n",
    "it doesn’t matter how the output was generated, as long as it is in the\n",
    "XML format.\n",
    "\n",
    "-   You can use Biopython to run BLAST over the internet, as described\n",
    "    in section \\[sec:running-www-blast\\].\n",
    "\n",
    "-   You can use Biopython to run BLAST locally, as described\n",
    "    in section \\[sec:running-local-blast\\].\n",
    "\n",
    "-   You can do the BLAST search yourself on the NCBI site through your\n",
    "    web browser, and then save the results. You need to choose XML as\n",
    "    the format in which to receive the results, and save the final BLAST\n",
    "    page you get (you know, the one with all of the\n",
    "    interesting results!) to a file.\n",
    "\n",
    "-   You can also run BLAST locally without using Biopython, and save the\n",
    "    output in a file. Again, you need to choose XML as the format in\n",
    "    which to receive the results.\n",
    "\n",
    "The important point is that you do not have to use Biopython scripts to\n",
    "fetch the data in order to be able to parse it. Doing things in one of\n",
    "these ways, you then need to get a handle to the results. In Python, a\n",
    "handle is just a nice general way of describing input to any info source\n",
    "so that the info can be retrieved using `read()` and `readline()`\n",
    "functions (see Section <span>sec:appendix-handles</span>).\n",
    "\n",
    "If you followed the code above for interacting with BLAST through a\n",
    "script, then you already have `result_handle`, the handle to the BLAST\n",
    "results. For example, using a GI number to do an online search:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIWWW\n",
    "result_handle = NCBIWWW.qblast(\"blastn\", \"nt\", \"8332116\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "If instead you ran BLAST some other way, and have the BLAST output (in\n",
    "XML format) in the file `my_blast.xml`, all you need to do is to open\n",
    "the file for reading:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "result_handle = open(\"my_blast.xml\", 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now that we’ve got a handle, we are ready to parse the output. The code\n",
    "to parse it is really quite small. If you expect a single BLAST result\n",
    "(i.e., you used a single query):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "blast_record = NCBIXML.read(result_handle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "or, if you have lots of results (i.e., multiple query sequences):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "blast_records = NCBIXML.parse(result_handle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Just like `Bio.SeqIO` and `Bio.AlignIO` (see\n",
    "Chapters \\[chapter:Bio.SeqIO\\] and \\[chapter:Bio.AlignIO\\]), we have a\n",
    "pair of input functions, `read` and `parse`, where `read` is for when\n",
    "you have exactly one object, and `parse` is an iterator for when you can\n",
    "have lots of objects – but instead of getting `SeqRecord` or\n",
    "`MultipleSeqAlignment` objects, we get BLAST record objects.\n",
    "\n",
    "To be able to handle the situation where the BLAST file may be huge,\n",
    "containing thousands of results, `NCBIXML.parse()` returns an iterator.\n",
    "In plain English, an iterator allows you to step through the BLAST\n",
    "output, retrieving BLAST records one by one for each BLAST search\n",
    "result:\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>TODO: should use an example with more than one record</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "result_handle = open(\"my_blast.xml\", 'r')\n",
    "blast_records = NCBIXML.parse(result_handle)\n",
    "blast_record = next(blast_records)\n",
    "print(blast_record.database_sequences)\n",
    "# # ... do something with blast_record"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Or, you can use a `for`-loop:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "result_handle = open(\"my_blast.xml\", 'r')\n",
    "blast_records = NCBIXML.parse(result_handle)\n",
    "for blast_record in blast_records:\n",
    "    print(blast_record.database_sequences)\n",
    "    # Do something with blast_record"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Note though that you can step through the BLAST records only once.\n",
    "Usually, from each BLAST record you would save the information that you\n",
    "are interested in. If you want to save all returned BLAST records, you\n",
    "can convert the iterator into a list:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "result_handle = open(\"my_blast.xml\", 'r')\n",
    "blast_records = NCBIXML.parse(result_handle)\n",
    "blast_records = list(blast_records)\n",
    "blast_records"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now you can access each BLAST record in the list with an index as usual.\n",
    "If your BLAST file is huge though, you may run into memory problems\n",
    "trying to save them all in a list.\n",
    "\n",
    "Usually, you’ll be running one BLAST search at a time. Then, all you\n",
    "need to do is to pick up the first (and only) BLAST record in\n",
    "`blast_records`:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "result_handle = open(\"my_blast.xml\", 'r')\n",
    "blast_records = NCBIXML.parse(result_handle)\n",
    "blast_record = next(blast_records)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "or more elegantly:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "result_handle = open(\"my_blast.xml\", 'r')\n",
    "blast_records = NCBIXML.parse(result_handle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "I guess by now you’re wondering what is in a BLAST record.\n",
    "\n",
    "The BLAST record class\n",
    "---------------------\n",
    "\n",
    "A BLAST Record contains everything you might ever want to extract from\n",
    "the BLAST output. Right now we’ll just show an example of how to get\n",
    "some info out of the BLAST report, but if you want something in\n",
    "particular that is not described here, look at the info on the record\n",
    "class in detail, and take a gander into the code or automatically\n",
    "generated documentation – the docstrings have lots of good info about\n",
    "what is stored in each piece of information.\n",
    "\n",
    "To continue with our example, let’s just print out some summary info\n",
    "about all hits in our blast report greater than a particular threshold.\n",
    "The following code does this:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "E_VALUE_THRESH = 0.04"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "result_handle = open(\"my_blast.xml\", 'r')\n",
    "blast_records = NCBIXML.parse(result_handle)\n",
    "\n",
    "for alignment in blast_record.alignments:\n",
    "    for hsp in alignment.hsps:\n",
    "        if hsp.expect < E_VALUE_THRESH:\n",
    "            print('****Alignment****')\n",
    "            print('sequence:', alignment.title)\n",
    "            print('length:', alignment.length)\n",
    "            print('e value:', hsp.expect)\n",
    "            print(hsp.query[0:75] + '...')\n",
    "            print(hsp.match[0:75] + '...')\n",
    "            print(hsp.sbjct[0:75] + '...')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Basically, you can do anything you want to with the info in the BLAST\n",
    "report once you have parsed it. This will, of course, depend on what you\n",
    "want to use it for, but hopefully this helps you get started on doing\n",
    "what you need to do!\n",
    "\n",
    "An important consideration for extracting information from a BLAST\n",
    "report is the type of objects that the information is stored in. In\n",
    "Biopython, the parsers return `Record` objects, either `Blast` or\n",
    "`PSIBlast` depending on what you are parsing. These objects are defined\n",
    "in `Bio.Blast.Record` and are quite complete.\n",
    "\n",
    "Here are my attempts at UML class diagrams for the `Blast` and\n",
    "`PSIBlast` record classes. If you are good at UML and see\n",
    "mistakes/improvements that can be made, please let me know. The Blast\n",
    "class diagram is shown in Figure\n",
    "\n",
    "![Class diagram for the Blast Record class representing all of the info\n",
    "in a BLAST report<span\n",
    "data-label=\"fig:blastrecord\"></span>](images/BlastRecord.png)\n",
    "\n",
    "The PSIBlast record object is similar, but has support for the rounds\n",
    "that are used in the iteration steps of PSIBlast. The class diagram for\n",
    "PSIBlast is shown in Figure.\n",
    "\n",
    "![Class diagram for the PSIBlast Record class.<span\n",
    "data-label=\"fig:psiblastrecord\"></span>](images/PSIBlastRecord.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>Deprecated BLAST parsers</font>\n",
    "-----------------------\n",
    "\n",
    "Bio.Blast.NCBIStandalone\n",
    "-----------------------\n",
    "The three functions for calling the \"legacy\" NCBI BLAST command line tools\n",
    "blastall, blastpgp and rpsblast were declared obsolete in Biopython Release\n",
    "1.53, deprecated in Release 1.61, and removed in Release 1.64. Please use\n",
    "the BLAST+ wrappers in Bio.Blast.Applications instead.\n",
    "\n",
    "The remainder of this module is a parser for the plain text BLAST output,\n",
    "which was declared obsolete in Release 1.54, and deprecated in Release 1.63.\n",
    "\n",
    "For some time now, both the NCBI and Biopython have encouraged people to\n",
    "parse the XML output instead, however Bio.SearchIO will initially attempt\n",
    "to support plain text BLAST output.\n",
    "\n",
    "Older versions of Biopython had parsers for BLAST output in plain text\n",
    "or HTML format. Over the years, we discovered that it is very hard to\n",
    "maintain these parsers in working order. Basically, any small change to\n",
    "the BLAST output in newly released BLAST versions tends to cause the\n",
    "plain text and HTML parsers to break. We therefore recommend parsing\n",
    "BLAST output in XML format, as described in\n",
    "section \\[sec:parsing-blast\\].\n",
    "\n",
    "Depending on which BLAST versions or programs you’re using, our plain\n",
    "text BLAST parser may or may not work. Use it at your own risk!\n",
    "\n",
    "### Parsing plain-text BLAST output\n",
    "\n",
    "The plain text BLAST parser is located in `Bio.Blast.NCBIStandalone`.\n",
    "\n",
    "As with the XML parser, we need to have a handle object that we can pass\n",
    "to the parser. The handle must implement the `readline()` method and do\n",
    "this properly. The common ways to get such a handle are to either use\n",
    "the provided `blastall` or `blastpgp` functions to run the local blast,\n",
    "or to run a local blast via the command line, and then do something like\n",
    "the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# result_handle = open(\"my_file_of_blast_output.txt\", 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Well, now that we’ve got a handle (which we’ll call `result_handle`), we\n",
    "are ready to parse it. This can be done with the following code:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# from Bio.Blast import NCBIStandalone\n",
    "# blast_parser = NCBIStandalone.BlastParser()\n",
    "# blast_record = blast_parser.parse(result_handle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This will parse the BLAST report into a Blast Record class (either a\n",
    "Blast or a PSIBlast record, depending on what you are parsing) so that\n",
    "you can extract the information from it. In our case, let’s just print\n",
    "out a quick summary of all of the alignments greater than some threshold\n",
    "value.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "E_VALUE_THRESH = 0.04\n",
    "for alignment in blast_record.alignments:\n",
    "    for hsp in alignment.hsps:\n",
    "        if hsp.expect < E_VALUE_THRESH:\n",
    "            print('****Alignment****')\n",
    "            print('sequence:', alignment.title)\n",
    "            print('length:', alignment.length)\n",
    "            print('e value:', hsp.expect)\n",
    "            print(hsp.query[0:75] + '...')\n",
    "            print(hsp.match[0:75] + '...')\n",
    "            print(hsp.sbjct[0:75] + '...')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "If you also read the section \\[sec:parsing-blast\\] on parsing BLAST XML\n",
    "output, you’ll notice that the above code is identical to what is found\n",
    "in that section. Once you parse something into a record class you can\n",
    "deal with it independent of the format of the original BLAST info you\n",
    "were parsing. Pretty snazzy!\n",
    "\n",
    "Sure, parsing one record is great, but I’ve got a BLAST file with tons\n",
    "of records – how can I parse them all? Well, fear not, the answer lies\n",
    "in the very next section.\n",
    "\n",
    "### Parsing a plain-text BLAST file full of BLAST runs\n",
    "\n",
    "We can do this using the blast iterator. To set up an iterator, we first\n",
    "set up a parser, to parse our blast reports in Blast Record objects:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# from Bio.Blast import NCBIStandalone\n",
    "# blast_parser = NCBIStandalone.BlastParser()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Then we will assume we have a handle to a bunch of blast records, which\n",
    "we’ll call `result_handle`. Getting a handle is described in full detail\n",
    "above in the blast parsing sections.\n",
    "\n",
    "Now that we’ve got a parser and a handle, we are ready to set up the\n",
    "iterator with the following command:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# blast_iterator = NCBIStandalone.Iterator(result_handle, blast_parser)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The second option, the parser, is optional. If we don’t supply a parser,\n",
    "then the iterator will just return the raw BLAST reports one at a time.\n",
    "\n",
    "Now that we’ve got an iterator, we start retrieving blast records\n",
    "(generated by our parser) using `next()`:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# blast_record = next(blast_iterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Each call to next will return a new record that we can deal with. Now we\n",
    "can iterate through these records and generate our old favorite, a nice\n",
    "little blast report:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# for blast_record in blast_iterator:\n",
    "#     E_VALUE_THRESH = 0.04\n",
    "#     for alignment in blast_record.alignments:\n",
    "#         for hsp in alignment.hsps:\n",
    "#             if hsp.expect < E_VALUE_THRESH:\n",
    "#                 print('****Alignment****')\n",
    "#                 print('sequence:', alignment.title)\n",
    "#                 print('length:', alignment.length)\n",
    "#                 print('e value:', hsp.expect)\n",
    "#                 if len(hsp.query) > 75:\n",
    "#                     dots = '...'\n",
    "#             else:\n",
    "#                 dots = ''\n",
    "#                 print(hsp.query[0:75] + dots)\n",
    "#                 print(hsp.match[0:75] + dots)\n",
    "#                 print(hsp.sbjct[0:75] + dots)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The iterator allows you to deal with huge blast records without any\n",
    "memory problems, since things are read in one at a time. I have parsed\n",
    "tremendously huge files without any problems using this.\n",
    "\n",
    "### Finding a bad record somewhere in a huge plain-text BLAST file\n",
    "\n",
    "One really ugly problem that happens to me is that I’ll be parsing a\n",
    "huge blast file for a while, and the parser will bomb out with a\n",
    "ValueError. This is a serious problem, since you can’t tell if the\n",
    "ValueError is due to a parser problem, or a problem with the BLAST. To\n",
    "make it even worse, you have no idea where the parse failed, so you\n",
    "can’t just ignore the error, since this could be ignoring an important\n",
    "data point.\n",
    "\n",
    "We used to have to make a little script to get around this problem, but\n",
    "the `Bio.Blast` module now includes a `BlastErrorParser` which really\n",
    "helps make this easier. The `BlastErrorParser` works very similar to the\n",
    "regular `BlastParser`, but it adds an extra layer of work by catching\n",
    "ValueErrors that are generated by the parser, and attempting to diagnose\n",
    "the errors.\n",
    "\n",
    "Let’s take a look at using this parser – first we define the file we are\n",
    "going to parse and the file to write the problem reports to:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import os\n",
    "# blast_file = os.path.join(os.getcwd(), \"blast_out\", \"big_blast.out\")\n",
    "# error_file = os.path.join(os.getcwd(), \"blast_out\", \"big_blast.problems\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now we want to get a `BlastErrorParser`:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# from Bio.Blast import NCBIStandalone\n",
    "# error_handle = open(error_file, \"w\")\n",
    "# blast_error_parser = NCBIStandalone.BlastErrorParser(error_handle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Notice that the parser take an optional argument of a handle. If a\n",
    "handle is passed, then the parser will write any blast records which\n",
    "generate a ValueError to this handle. Otherwise, these records will not\n",
    "be recorded.\n",
    "\n",
    "Now we can use the `BlastErrorParser` just like a regular blast parser.\n",
    "Specifically, we might want to make an iterator that goes through our\n",
    "blast records one at a time and parses them with the error parser:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# result_handle = open(blast_file)\n",
    "# iterator = NCBIStandalone.Iterator(result_handle, blast_error_parser)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We can read these records one a time, but now we can catch and deal with\n",
    "errors that are due to problems with Blast (and not with the parser\n",
    "itself):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# try:\n",
    "# next_record = next(iterator)\n",
    "# except NCBIStandalone.LowQualityBlastError as info:\n",
    "# print(\"LowQualityBlastError detected in id %s\" % info[1])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
